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Abstract 

It  is  very  difficult  to  standardize  programs  in  mathematical 
physics. 

^received  on  August  id,  T981 

1)  In  addition,  comrades  Liu  Zhiping  ahd  Li  Yingtao  of  the 
computing  center,  Gu  Peizi,  Wu  Zhiyuan,  Xiao  Zuoshi  and  Xu 
Feng  of  the  third  chemical  engineering  institute,  and  Zhou 
I  Xiaowei  of  Shanxi  Qinghua  Electrical  Machine  Institute  also 
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,  participated  in  this  work. 
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However,  it  is  not  poaslble  to  compile  an  application  software  or  /130 
program  for  a  specific  type  of  problem.  We  prim^arily  referred  to 
the  "HELP  Code"^^^'^to  compile  a  LTDL  program^ ^  in  BCY  language. 

It  is  capable  of  solving  plasticity  and  elasticity  problems  of  a 
two  dimensional  nonstationary  fluid  (planar  or  axlsymmetry 
problems),  including  multi-materials  (metals  and  non-metals).  A 
great  deal  of  calculation  was  performed  to  solve  problems  such  as 
metal  projectile  penetration  into  steel  target,  high  speed 
collision,  surface  explosion  of  high  energy  explosive,  spherical 
explosion,  collision  of  two  rods,  etc.  Very  useful  results  were 
obtained. 

This  paper  describe^he  capabilities  and  focal  points 

of  this  software. In  the  first  section,  we  will  introduce  the 
fundamentals  of  thus  software,  i.e.  program  modulation  based  on 


^ 
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split  operator  method.  In  the  second  section,  the  key  of 
software  design  -  calculation  of  mixed  cell  flow  -  is  described. 

The  third  section  describes  the  determination  of  the  problem, 
i.e.  setting  up  the  region  and  boundary.  The  fourth  section 
discusses  techniques  to  divide  meshes.  Energy  verification  is 
Introduced  in  section  five.  Several  examples  and  figures  are 
shown  in  section  six. 

§1.  Split  Operator  Method  -  Modulation  of  Numerical  Calculation 

In  order  to  compile  a  software  for  a  specific  type  of 
physical  problem  which  is  affected  by  many  factors,  the  split 
operator  method  may  be  used  in  approximation.  The  effect  of 
several  factors  is  split  into  the  superposition  of  effects  of 
several  single  factors.  A  numerical  calculation  program  based  on  the 
split  factor  method  is  modular.  A  new  physical  problem  may  be 
calculated  by  suitably  selecting  the  relevant  operators 
(corresponding  to  several  modules)  and  to  match  them  with 
boundary  conditions  and  region  of  interest. 

The  split  operator  method  was  initially  used  to  lower  the 
number  of  dimensions.  For  example,  it  can  reduce  a  three- 
dimensional  problem  to  three  one-dimensional  problems,  i.e. 

m,  —  L(*,,  •  L,,  •  L,,  •  «. 

Moreover,  the  split  operator-method  can  be  applied  to  different 
physical  effects.  For  instance,  the  N-S  equation  can  be  considered 
as  the  superposition  of  the  invlscous  operator  and  viscous 
operator  ,  i.e. 

■«  ™  •  L,  • 

Similarly,  the  plasticity  and  elasticity  operator  of  a  fluid  may 
be  split  into  the  pressure  effect  operator  ,  stress  effect 
operator  L^,  and  the  transport  effect  operator  i.e. 

■i  •  A,  •  L,  •  L,  • 

Each  operator  is  separated  to  become  independent  modules  when  the 
numerical  calculation  program  is  compiled.  For  instance,  the 
LTDL  program  was  compiled  based  on  this  principle.  This 


r 


program  is  very  flexible.  For  example,  a  fluid  can  be  calculated 


by  shutting  off  program  modules  corresponding  to  in  LTDL.  If 


there  is  only  one  fluid  material  and  it  includes  a  fl?ee  boundary, 

r2 

i.e.  the  FLIC  method  *  ,  then  all  modules  corresponding  to  L, , 
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L2  and  are  used,  which  is  the  so-called  HELP'-  method. 


§2.  Key  Of  Software  Design  -  Mixed  Cell  Flow  Calculation 


energy. 


The  most  complex  portion  of  the  software  to  calculate 
plasticity  and  elasticity  of  two  dimensional  nonstationary  fluids 
is  mixed  cell  flow  calculation  (cells  at  the  interface  of 
materials).  The  key  is  to  determine  the  rules  to  calculate  the  /131 

area  of  material  Interface  which  cuts  the  Euler  mesh.  Obviously, 
the  direction  of  flow,  the  shape  of  interface,  the  number  of 
interfaces  in  a  section  cell  can  vary  widely.  To  use  a  set  of 
rules  to  calculate  the  area  where  the  Interface  cuts  into  a  cell, 
the  density ‘of  the  transport  mesh  and  the  transport  rate  at  the 
interface  will  enable  us  to  calculate  the  mass,  momentum  and 
energy  transport  between  a  mixed  cell  and  its  neighbors.  Then, 
the  mechanical  parameters  of  a  mixed  cell  at  any  time  can  be 
calculated  by  using  laws  of  conservation  of  mass,  momentum  and 
energy.  The  logic  is  very  complicated  and  includes  the  following  steps: 

1.  First,  the  interface  of  materials  is  defined. 

Initially,  a  set  of  massless  tracing  points  (their  coordinates 
are  expressed  by  the  lattice  coordinates).  These  points  are 
labelled  so  that  the  material  is  always  on  the  left  side  of  any 
two  neighboring  points.  The  points  at  the  interface  coincide  and 
move  in  opposite  directions.  The  position  of  a  moving  point  can 
be  determined  bjf  a  weighted  mesh  speed  area  calculation.  The 
broken  line  between  two  neighboring  points  is  the  interface. 

2.  The  area  under  which  the  interface  Intersects  n  sides  of 
a  lattice  is  calculated  according  to  the  following  three 
situations : 

1)  The  interface  only  Intersects  a  side  of  the  lattice  once. 

In  Figure  1,  the  interface  EBCDF  between  materials  (1)  and 

(2)  Intersects  with  a  side  of  lattice  K  once.  The  area  is 
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calculated  by  linear  interpolation.  The  Intersect  E  between,  the 
Interface  and  the  top  side  of  cell  K  is  calculated  from  the 
coordinates  of  two  neighboring  points  on  the  top  edge  of  cell 
K.  Then,  based  on  the  rule  that  the  material  is  always  on  the  left 
of  the  line  connecting  two  neighboring  points,  we  can  then  determine 
that  the  area  under  which  material  (1)  Intersects  with  the  top 
side  of  cell  K  is  the  convolution  of  AE.  The  area  under  which 
material  (2)  Intersects  with  the  top  side  of  cell  K  is  the 
convolution  of  EP. 

11)  The  Interface  Intersects  the  lattice  several  times. 

In  Figure  1,  the  Interface  of  materials  (1)  and  (2) 
intersects  with  the  right  side  of  lattice  K  three  times.  The 
area  of  intersection  is  obtained  by  adding  the  area  of  each 
intersection  and  by  comparing  that  area  with  the  entire  area  of 
cell  K,  Specifically,  method  i)  is  used  along  the  direction  of 
the  tracing  points  of  the  interface  of  material  (1)  to  first 
calculate  the  area  BA  and  then  the  second  area  CH.  Furthermore, 

CH  is  added  to  BA.  If  BA  +  CH  is  less  than  AH,  then  the  sum 
remains  unchanged.  Finally,  the  area  of  the  third  intersection 
DA  is  calculated  and  added  to  BA  and  CH.  If  BA  +  CH  +  DH  >  AH, 
then  AH  is  subtracted  from  the  sum  to  obtain  BA  +  DC,  which  is 
the  area  of  three  intersections  between  Interface  (1)  and  the 
right  side  of  cell  K.  The  derivation  for  n  intersections  can  be 
obtained  in  analogy. 


Figure  1 


Ill)  There  Is  no  Intersection  between  the  interface  and  n 
sides  of  the  lattice 

In  Figure  1,  the  Interface  of  materials  (1)  and  (2)  does,  not 
Intersect  the  bottom  and  left  sides  of  cell  K.  From  the  figure 
we  know  that  the  area  of  Intersection  between  the  Interface  of 
material  (1)  and  the  left  and  bottom  sides  of  cell  K  Is  zero. 

That  of  the  Interface  of  material  (2)  Is  the  entire  cell.  This 
Is  achieved  by  "pre-set"  (Interface  entering  cell  K)  and  "re-set" 
(Interface  leaving  cell  K)  In  the  program.  Before  we  describe 
this  method,  let  us  first  Introduce  the  four  directions 
connecting  the  tracing  points  on  the  top  and  right  sides  of 
lattice  K  (See  Figure  2). 


Figure  2 

(1)  entering  cell  K,  leaving  cell  KR. 

(2)  leaving  cell  K,  entering  cell  KR. 

(3)  entering  cell  K,  leaving  cell  KA. 

(4)  leaving  cell  K,  entering  cell  KA. 

"Pre-set"  rules.  When  the  Interface  of  a  material 

Intersects  with  cell  K  (on  the  top  or  right  side)  n  times  (nsl), 
If  the  area  of  Intersection  with  cell  K  Is  less  than  the  total 
area  of  the  whole  cell,  other  sides  of  cell  K  will  not  be  "pre¬ 
set".  Otherwise,  other  sides  of  the  cell  are  pre-set  clockwise 
starting  from  the  Intersect  until  reaching  the  side  of 
Intersection  again. 


"Re-set"  rules.  When  a  material  Interface  Intersects  n 
times  (n^1)  with  a  side  of  cell  K  (top  or  right)  and  leaves  cell 
K,  If  the  cumulated  area  of  Intersection  Is  less  than  the  total 
area,  then  other  sides  of  cell  K  will  not  be  re-set.  Otherwise, 

It  will  be  re-set  clockwise  from  the  Intersect.  The  partial  area 
of  other  side  Is  zero  until  the  Interface  Intersects  with  the  lattice 
side. 

Let  us  use  Figure  1  to  explain  this  situation.  If  we  scan 
along  the  direction  of  tracing  points  for  material  (1),  when  the 
line  connecting  two  such  points  for  material  (1)  Intersects  with 
the  top  of  cell  K  once,  then  the  right,  bottom  and  left  sides  of 
cell  K  are  pre-set  to  the  total  area  according  to  the  pre-set 
rule.  When  the  Interface  of  material  (1)  Intersects  with  the 
right  side  of  cell  K  (at  B)  and  leaves  cell  K,  the  left  side  Is 
re-set  to  zero  area.  Because  the  top  side  of  cell  K  already 
Intersects  with  the  Interface  of  material  (1),  therefore.  It  Is 
not  re-set  again.  When  the  Interface  of  material  (1)  Intersects 
with  the  right  side  of  cell  K  three  times  (at  points  B  and  C) 
before  entering  cell  K,  because  the  cumulative  area  Is  smaller 
than  the  total  area  to  the  right  of  cell  K,  then  the  bottom  and 

left  sides  are  "re-set"  to  the  total  area.  When  the  Interface  of 

material  (1)  Intersects  with  the  right  side  of  cell  K  three  times 
(at  points  B,  C  and  D),  because  the  cumulative  area  Is  larger 
than  the  total  area  to  the  right  of  cell  K,  then  the  bottom  and 

left  sides  are  "re-set"  to  zero.  This  appears  to  be  correct  from 

Figure  1.  Similar  "pre-set"  and  "re-set"  rules  may  be  used  to 
determine  that  the  area  of  Intersection  between  the  Interface  of 
material  (2)  and  the  bottom  and  left  sides  of  cell  K  Is  the  total 
area. 

In  summary,  based  on  the  direction  of  the  tracing  points  of 
the  "pre-set"  and  "re-set"  rules,  we  can  calculate  any  partial 
area  of  Intersection  between  an  Interface  and  a  side  of  the 
lattice . 


§3.  Determination  of  the  Problem ,  Region  and  Boundary  Condition 


The  LTDL  program  Is  capable  of  calculating  mutual 
Interaction  among  various  materials.  A  material  may  also  Include 
other  types  of  materials  (distributed  In  different  regions).  The 
shape  of  the  material  may  be  spherical,  hemispherical, 
cylindrical  and  conical,  as  well  as  their  combinations.  However, 
they  must  be  axlsymmetrlc  or  planar.  The  boundary  conditions 
Include  symmetric  axis,  rigid  boundary,  transport  boundary,  free 
boundary  and  Interfaclal  boundary.  With  the  exception  of 
transport  surface  and  symmetric  axis,  other  Interfaces  have 
massless  tracing  points. 

Let  us  assume  that  the  area  of  calculation  Is  the 
rectangular  region  ABCD  on  the  (r,  Z)  surface  In  a  cylindrical 
coordinate,  as  shown  in  Figure  3.  Initially,  the  boundaries  of 
material  blocks  (1),  (2)  and  (3)  are  marked  by  massless  points. 
Thus,  the  lattice  of  calculation  may  ha  a  pure  cell  of  a  single 
material  or  may  contain  two  or  three  mixed  materials.  Initial 
values  of  eight  mechanical  parameters  such  as  u  (velocity 
component  in  r-direction) ,  v  (velocity  component  In  Z-direction) , 
e  (specific  Internal  energy),  p  (density),  (stress 

components)  and  p  (pressure)  are  given  as  the  Initial  condition. 
Thus,  based  on  the  problem  (Including  the  Initial  condition  of 
the  material),  shape,  boundary  condition,  equation  of  state  and 
strength  formula,  the  software  can  be  rationally  used  to  perform 
the  corresponding  calculation.  In  addition  to  the  eight 
parameters  and  the  dynamic  boundary  of  the  material,  It  Is  also 
necessary  to  provide  strength,  temperature,  positive  stress  and 
strain  rate. 


The  determination  of  area  of  calculation  and  boundary  type 
can  be- explained  with  the  following  example. 

1.  Axial  Collision  of  Two  Semi-infinite  Metallic  Bodies 


Figure  4  Figure  5  Figure  6 

In  reality,  it  is  an  axial  one-dimensional  problem  and  the 
area  of  calculation  is  ABCD  (Figure  4).  AB,  BF,  FC  and  CD  are 
transport  boundaries.  AED  is  the  axis  of  symmetry.  There  is  no 
need  for  tracing  points.  EF  is  the  interface  between  blocks  (1) 
and  (2)  which  requires  tracing  points. 


2.  Axial  Collision  of  Two  Finite  Thickness  Infinite  Width 
Materials 

This  is  also  a  one -dimensional  axial  problem.  The  area  of 
calculation  is  GHJI  (Figure  5).  AB  and  CD  are  free  boundaries. 

EF  is  an  interface.  They  all  need  tracing  points.  BF  and  FC  are 
transport  boundaries  and  GAEDl  is  the  axis  of  symmetry. 

3.  Axial  Collision  of  Two  Finite  Thickness  and  Width 
Materials 

The  area  of  calculation  of  this  two-dimensional  axlsymmetric 
problem  is  GHJI  (Figure  6).  AB,  BF,  FC  and  CD  are  free 
boundaries.  EF  is  an  interface.  They  should  have  tracing 
points.  AED  is  the  axis  of  symmetry. 

4.  High  Speed  Collision  of  Copper  Sphere  with  a  Semi- 
infinite  Body  (Figure  7) 

Because  of  symmetry,  it  is  only  necessary  to  calculate  half 
of  the  area  (Figure  8).  Because  the  target  is  a  semi-infinite 
body.  Therefore,  CD  and  DB  are  transport  boundaries .AB  and  AHG 
are  free  boundaries.  CAGE  is  the  axis  of  symmetry.  The  area  of 
calculation  is  CDFE. 
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5.  Jet  Penetrating  Steel  Target 

If  the  target  has  a  finite  thickness,  then  the  back  surface 
of  the  target  CD  is  a  free  boundary.  CD  is  a  transport  boundary 
if  the  target  is  infinitely  thick.  GB  is  a  free  boundary.  Let 
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US  assume  that  the  jet  is  Infinitely  long,  then  EH  is  a  transport 
boundary  and  GH  is  a  free  boundary.  AG  is  an  interface  and  IDAE 
is  the  axis  of  symmetry.  When  DC  and  CB  are  free  surfaces,  the 
area  of  calculation  is  IJKE.  When  DC  and  CB  are  transport 
surfaces,  the  area  of  calculation  is  DCFE. 

6.  Surface  Explosion  of  High  Energy  Explosive 

Let  us  assume  that  there  is  a  high  energy  explosive  at  GFED 
on  the  surface  of  a  solid  body.  ABCEFG  is  a  solid.  Outside  is 
the  atmosphere.  Then,  it  is  necessary  to  make  AB  and  BC  be 
transport  surfaces,  CD  be  a  free  surface,  GFED  be  an  interface 
and  AGDH  be  the  axis  of  symmetry. 
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Figure  9 


Figure  10 

1.  solid 

2.  air 


§4.  Expansion  of  Area  of  Calculation  -  Techniques  for  Lattice 
Reconf igur at ion 


The  special  feature  of  a  nonstationary  motion  is  that  the 
pattern  of  motion  varies  with  time.  However,  due  to  limited 


computer  storage  space,  the  number  of  Euler  lattices  cannot  be 
too  many.  After  a  certain  time,  the  physical  picture  has  filled 
the  entire  area  of  calculation.  Then,  lattices  can  be  combined 
by  reconfiguration  to  expand  the  space  of  calculation.  Lattice 
combination  can  be  done  In  one  or  two  directions.  For  example, 
the  program  Is  calculated  to  some  extent,  perturbation  may  have 
already  affected  all  the  cells.  Thus,  the  original  100  x  50 
cells  are  consolidated  Into  50  x  25  or  100  x  25  cells.  Thus,  3/4 
or  1/2  of  the  storage  space  can  be  liberated  to  extend  the  area 
of  calculation. 

The  mechanical  parameters,  positions  of  tracing  points  at 
the  Interface,  positions  of  lattice  lines  and  cell  sizes  must  be 
redetermined  after  reconfiguration.  The  calculation  can  be 
performed  In  four  steps: 

1.  Calculate  the  position  of  each  lattice  line,  the 
coordinate  of  each  tracing  point  and  the  axial  and  radial 
dimensions  of  each  lattice. 

2.  Calculate  mechanical  parameters  In  reconfigured  cells. 

Let  us  assume  that  cells  L  and  ^  are  combined  Into  cell  K.  The 
principle  of  recombination  Is  as  follows:  the  mass,  axial  and 
radial  velocity,  and  specific  Internal  energy  In  cell  K  are 
calculated  from  the  conservation  of  mass,  momentum  and  energy  In 
cells  L  and  M.  The  stress  derivation  of  cell  K  Is  calculated  by 
weighted  mass  average  of  cells  L  and  M.  If  the  combined  cell  Is 
a  mixed  cell,  then  the  mass,  velocity  and  specific  Internal 
energy  of  the  cell  are  obtained  by  laws  of  conservation.  The 
density  of  each  material  In  cell  K  can  be  obtained  by  weighted 
density  mass  of  the  same  method  In  cells  L  and  M. 

3.  New  Interface  and  center  tracing  points  are  added  to  the 
new  lattice  In  order  to  calculate  the  position  of  the  lattice  line 
and  the  size  of  the  cell.  In  addition,  mechanical  parameters  In 
the  new  cells  are  defined  by  those  In  the  boundary  cells. 

4.  Based  on  the  mixed  cell  pressure  Iteration  principle 
(see  reference  Cl]  for  details)  and  equation  of  state,  the 
pressure  of  each  cell  (pure  and  mixed)  Is  calculated.  Finally,  /135 
the  time  Interval  Is  recalculated. 


In  regard  to  the  treatment  of  a  free  boundary,  the  technique 
is  shown  in  detail  In  reference  [53. 

§5.  Energy  Verification 

A  stationary  problem  Is  often  associated  with  Iteration. 
Therefore,  some  local  error  will  not  affect  the  entire  picture. 

At  best,  a  few  more  Iterations  are  needed.  However,  the  problem 
of  nonstationary  calculation  Is  different.  It  Is  necessary  to 
check  for  energy  conservation  at  each  level  (each  time  Interval) 
for  every  cell.  Otherwise,  the  result  cannot  be  trusted. 

Initially,  the  total  energy  of  the  system  ETH  Is  calculated 
according  to  the  lattice  number.  The  theoretical  energy  ETH*  of 
the  system  at  the  nth  level  Is  obtained  by  adding  to  or 
subtracting  from  ETH  with  the  energy  Input  (work  done  by  the 
outside)  or  energy  output  (work  done  by  the  system)  across  the 
transport  boundary  In  a  time  Interval  At.  The  total  energy  at 
the  nth  level  Is  ESUM.  It  Is  required  that  the  relative  error 
should  be  within  the  range  DMIN,  l.e. 

[ETH*  -  ESUM 
I  ETH* 

This  formula  must  be  satisfied  to  demonstrate  that  the 
calculation  Is  correct.  Otherwise,  It  has  problems  and  the 
calculation  Is  terminated. 

If  this  energy  verification  Is  used  for  each  cell,  the 
requirement  Is  even  more  rigorous. 

§6.  Examples  of  Numerical  Calculation 

1 .  Jet  Penetrating  Steel  Target  Plate 
Initial  conditions  are:  target  radius  17.44mm,  thickness 
infinite  (actual  thickness  is  18.49mm),  p^  «  7.85g/cm*  ,  p  «  latm, 
e^O. 6369426  x  10*ergs/g,  u  *  v  *  0.  Jet  (copper)  radius  is  1mm. 
The  Jet  Is  Infinitely  long  (8mm  In  calculation),  «  8.60g/cm’ , 
p  ■  latm,  e*.  *  0.5813980  x  10*ergs/g,  u  »  0,  v  ■  v*.  ■  850m/s. 


Figure  11.  Jet  Penetration  of  Steel  Plate  Ended  at  t  «  2  s  and 
Stationary  Penetration  Began' 

1.  jet 

2.  steel 

2.  Coppy  Ball  Colliding  with  Steel  Target  at  2kin/sec. 
Figures  8,  12  and  13  show  the  external  shape  at  t  >  0, 

0.12247  and  0.4487  microseconds, 

3.  Spherical  Explosion  in  Vacuum^^^ 

Figures  14  and  15  show  the  external  profile  at  t  >  0  and  0.3 
usee,  respectively. 


Figure  12 

1.  n  =  40 

2.  t  =  0.12247)tsec 
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Figure  14 


Figure  13 

1 .  n  “  121 

2.  t  »  0.4487Msec 


Figure  15 
1 .  T  «  0.3usec 


Figure  16 
1 .  T  *  70ttsec 
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Figure  18 

1 .  explosive 

2.  charge  mask 


.Figure  17 
1 .  T  *  0.4601  tisec 


Figure  19.  Shapes  and  Posl> 
tlons  of  Explosion 
Product  and  Charge 
Mask  at  t  s6.468)isec 

1.  t  >  6.468)isec 


Figure  20.  Penetration  of  Shield  In  5/«ec 

4.  Surface  Explosion  of  High  Energy  Explosive 

Figure  10  shows  the  profile  at  t  »  0.  Figure  16  Is  the 
pattern  at  70ms.  The  expansion  of  gas  explosive  products  can  be 
seen  In  Figure  16. 

5.  Collision  of  Two  Rods 

Figures  6  and  17  are  the  profiles  at  t  »  0  and  0.4601ms, 
respectively. 

6.  Calculation  of  Large  Conical  Angle  Tumbling  Projectile 

Figures  18  and  19  are  the  profiles  at  t  >  0  and  6.486ms, 
respectively. 

7.  Calculation  of  Cylindrical  Steel  Projectile  Hitting 
Steel  Target  at  900m/ sT^-^ 

Figure  20  Is  the  profile  at  t  »  5ms. 
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